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Many  background  error  correlation  (BEC)  models  in  data  assimilation  are  for¬ 
mulated  in  terms  of  a  positive-definite  smoothing  operator  B  which  simulates 
the  action  of  correlation  matrix  on  a  vector  in  state  space.  To  estimate  the  ef¬ 
ficiency  of  such  approach,  numerical  experiments  with  the  Gaussian  and  spline 
models 


have  been  conducted.  Here  I  is  the  identity  operator  and  u  is  the  diffusion 
tensor,  whose  spatial  variability  is  derived  from  the  forecast  field  and  m  is  the 
spline  approximation  order. 

Performance  of  these  BEC  representations  are  compared  in  the  frame¬ 
work  of  numerical  experiments  with  real  3dVar  data  assimilation  into  the  Navy 
Coastal  Ocean  model  (NCOM)  in  the  Western  Tropical  Pacific.  It  is  shown  that 
both  BEC  models  have  similar  forecast  skills  over  a  two-month  time  period, 
whereas  the  second-order  spline  model  is  several  times  more  efficient  compu¬ 
tationally  if  the  cast  function  is  minimized  in  the  state  space. 
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1.  Introduction 

In  recent  years,  heuristic  BEC  modelling  has  become  an  area  of  active 
research  in  data  assimilation  (DA).  This  interest  has  been  fueled  by  de¬ 
velopment  of  the  ensemble  DA  techniques  and  rapid  increase  of  the  data 
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streams,  driven  by  remotely  sensed  observations  from  satellites  and  intro¬ 
duction  of  new  autonomous  observational  platforms  in  oceanography.  Tra¬ 
ditional  BEC  models  based  on  the  explicit  definition  of  the  error  covariances 
by  the  families  of  parameterized  correlation  functions1,2  tend  to  loose  com¬ 
putational  efficiency  with  the  growth  of  the  number  of  observations  and 
with  the  necessity  to  introduce  more  complex  observation  operators  into 
the  DA  algorithms.  Because  of  that,  there  is  a  growing  tendency  to  esti¬ 
mate  error  correlations  directly  in  the  model  space  (MS)  or  in  its  subspaces 
spanned  by  the  appropriately  selected  basis  functions.3  5 . 

Of  particular  interest  are  the  MS  correlation  models  based  on  positive 
functions  of  the  diffusion  operator  D  =  Vi/V,  such  as  the  exponent  or 
the  inverse  of  its  binomial  approximation6,7.  Both  types  of  models  were 
extensively  used  in  many  applications8'11  because  of  their  convenience  and 
ease  of  numerical  implementation.  Another  advantage  is  their  flexibility  in 
approximation  of  inhomogeneous  and  anisotropic  covariance  functions12,13. 
Numerically,  these  models  are  implemented  by  integrating  the  diffusion 
equation  (DE)  using  either  explicit  or  implicit  scheme. 

The  computational  cost  of  the  DE  integration  by  explicit  methods  in¬ 
creases  substantially  when  the  local  decorrelation  scale,  pc ,  becomes  larger 
than  the  model  grid  step  6x.  This  is  because  the  minimum  number  of  mul¬ 
tiplications  by  D,  is  proportional  to  (pc/6x)2  -  a  constraint  imposed  by 
numerical  stability  of  the  integration.  In  such  situations  it  may  be  advanta¬ 
geous  to  employ  implicit  integration  schemes11,14,  which  tend  to  converge 
fast  enough  to  deliver  considerable  computational  gain.  Additional  gain  can 
be  obtained  if  the  implicit  approximation  is  implemented  within  the  model 
space  (MS)  formulation.  In  this  case,  the  iterative  solution  of  the  system  in 
data  space  (DS)  which  embeds  an  iterative  cycle  of  the  implicit  scheme,  is 
no  longer  needed. 

This  study  compares  the  forecast  skill  and  computational  cost  of  two 
BEC  models:  The  first  model  (Coo)  is  described  by  the  propagator  of  the  DE 
and  implemented  numerically  by  its  explicit  integration;  the  second  BEC 
model  (Cm)  is  defined  by  the  inverse  of  a  mth-order  binomial  of  D,  that  ap¬ 
proximates  Coo  and  can  be  interpreted  as  a  result  of  m-step  DE  integration 
with  the  implicit  scheme.  Assimilation  experiments  were  performed  using 
both  DS  and  MS  formulations  of  Cm  and  a  realistic  regional  ocean  model 
with  real  data.  It  is  shown  that  in  certain  situations  it  is  computationally 
advantageous  to  employ  the  second  (spline)  model  combined  with  the  MS 
solution  to  the  normal  equation. 
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2.  Covariance  modelling  in  MS 

2.1.  MS  and  DS  approaches  in  SdVar 

In  its  basic  formulation,  the  3dVar  analysis  determines  the  optimal  model 
state  increment  x  that  minimizes  the  following  cost  function, 

J(x)  =  l-  [xtB-1x  +  (Hx  —  d)TR  "'(Hx  —  d)]  -4  min.  (1) 

Here  B  is  the  BE  covariance  matrix  and  d  is  the  innovation  vector 

d  =  y  —  Hxb 

where  y  denotes  observations,  x*,  is  the  background  model  state,  H  is  the 
observation  operator  (linearized)  in  the  vicinity  of  x&,  and  R  is  the  covari¬ 
ance  matrix  of  observation  errors.  To  simplify  the  notation,  the  variables  in 
both  model  and  data  spaces  are  non-dimensionalized  by  x  4—  2x  and 
d  «—  R_1/  d;  where  C  is  the  (diagonal)  background  error  variance  matrix. 
In  order  to  keep  J  invariant,  the  matrices  B,  and  H  are  non-dimensionalized 
by  B  C~1/2BC"1/2;  H  f- R_1/2HC1/2. 

The  cost  function  (1)  is  minimized  by  solving  the  normal  equation  which 
sets  the  gradient  of  J  equal  to  zero: 

(B'  +  HTH)x  =  HTd  (2) 

so  that  the  solution  to  the  normal  equation  is: 

x  =  (B  1  +HTH)-1HTd  (3) 

Solving  equation  (2)  for  the  model  state  increment  x  is  the  basic  tool  of 
3dVar  analysis.  If  B_1  has  full  rank,  the  solution  (3)  is  unique  and  can  be 
rewritten  in  the  dual  form15: 

x  =  BHT(HBHT  +  I)"1d  (4) 

which  is  often  called  the  DS  solution  to  the  variational  problem  (1).  Note 
that  if  B_1  does  not  have  full  rank,  defining  B  as  its  generalized  inverse 
does  not  guarantee  that  solution  (4)  will  coincide  with  the  solution  (3)  of 
the  original  minimization  problem.  This  is  because  the  DS  solution  (4)  is 
always  orthogonal  to  the  null  space  of  B,  whereas  in  general,  the  minimizer 
(2)  of  (1)  is  not  constrained  by  this  condition. 

It  should  also  be  noted  that  the  majority  of  the  BEC  models  are  based 
on  direct  computation  of  the  matrix  elements  of  B  from  experimental  data, 
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and  therefore,  require  3dVar  solutions  in  the  forms,  involving  B  These  are 
the  DS  solution  (4),  or  the  MS  solution,  preconditioned  by  B  =  VB  6’17: 

x  =  B(l  +  BIHTHB)-|BHTd  (5) 

Therefore,  the  basic  MS  solution  (3)  has  been  used  in  practice  more  rarely 

for  two  reasons:  a)  it  requires  solving  the  linear  system  in  MS,  which  has 
many  more  dimensions  than  the  DS;  and  b)  estimation  of  B  from  the  data 
is  more  straightforward  than  estimation  of  B”1. 

2.2.  The  Gaussian  and  spline  BEC  models 

The  two  types  of  BEC  models  considered  here  are  based  on  the  polynomials 
of  D.  The  major  idea  is  to  model  the  resulting  action  of  the  BEC  operator 
B  on  a  vector  x  by  integrating  the  corresponding  diffusion  equation 

dx  ^  1„  „ 

dt  =  °X  S  2Vl/Vx  (6) 

for  a  certain  ’’time  period”  T,  thus  setting  B  =  expTD. 

The  diffusion  tensor  is  is  represented  by  3x3  positive-definite  matrices 
whose  entries  depend  on  the  coordinates  x  in  physical  space.  The  eigen¬ 
values  A^,i  =  1,...,3  of  isT  are  all  positive,  have  the  dimension  of  length 
squared,  and  in  the  homogeneous  case  (is  —  const)  they  are  naturally  inter¬ 
preted  as  the  squares  of  the  decorrelation  scales  pt  in  the  directions  of  the 
respective  eigenvectors  of  is.  In  the  inhomogeneous  case,  the  decorrelation 
scales  are  defined  locally  in  a  similar  manner,  whereas  the  integration  time 
T  plays  the  role  of  a  global  scaling  parameter  for  the  distribution  of  p2(x). 
Therefore,  setting  the  value  of  T  is  equivalent  to  specifying  the  square  of 
the  mean  decorrelation  scale  p  for  a  given  distribution  of  is(x).  Throughout 
the  remainder  of  this  paper,  we  keep  in  mind  this  equivalence  and  replace 
T  with  p2  where  appropriate. 

Numerically,  the  action  of  the  Gaussian  BEC  operator  exp(TD)  is  usu¬ 
ally  represented  by  integrating  (6)  with  an  explicit  time-stepping  scheme, 
xt+<5t  __  xe  +  (^Dx*,  such  that  the  result  of  multiplication  of  a  vector  x°  by 
B  is 


xT  =  Bx°  = 


TD 

1  + - 

n 


x°  «  exp[p2D]x° 


(7) 


where  n  =  T/St  is  the  total  number  of  time  steps.  Expression  (7)  shows  that 
numerically,  the  Gaussian  BEC  model  is  a  high-order  polynomial  in  D.  In 
data  assimilation  problems  the  n-step  ’’time  integration”  (7)  is  embedded 
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in  the  iterative  loop  that  solves  linear  equations,  whose  system  matrices  are 
under  the  inversion  signs  in  either  (4)  or  (5).  Therefore,  reducing  n  (increas¬ 
ing  St)  provides  major  computational  savings.  The  minimum  value  for  n  is 
limited,  however,  by  the  stability  condition  which  constrains  eigenvalues  of 
the  operator  I  +  p2D/n  in  (7)  not  to  exceed  1  in  magnitude: 

n  >  ip2  A  (8) 

where  A  is  the  absolute  value  of  the  largest  eigenvalue  of  D.  Numerically, 
the  minimum  value  of  n  in  3d  is  proportional  to  the  square  of  the  largest 
ratio  p  between  decorrellation  scale  and  the  local  grid  step  taken  over  the 
entire  grid.  In  realistic  applications,  p  may  easily  exceed  10,  substantially 
increasing  the  cost  of  computing  the  action  of  B  on  a  vector. 

For  p  >  10  the  computational  burden  can  be  reduced  by  considering  a 
spline  BEC  model 


B 


m 


m 


exp[p2D], 


(9) 


which  specifies  the  inverse  BEC  as  a  polynomial  in  the  powers  of  —  D  and 
converges  to  the  Gaussian  model  as  m  — >  oo  (Fig.  1). 


Fig.  1.  Normalized  spectra  of  the  Gaussian  (m  =  oo)  and  spline  BEC  operators  with 
different  approximation  order  m  (a)  and  the  corresponding  correlation  functions  (b). 
Horizontal  axis  is  normalized  by  the  correlation  radius.  The  dashed  line  shows  correlation 
function  used  in  the  experiments  with  NCODA  Cd  model  (see  Section  3.2). 


The  BEC  operator  in  (9)  can  be  implemented  numerically  in  two  ways, 
distinguished  by  the  order  of  the  operations  of  inversion  and  raising  to 
the  rnth  power.  The  first  method  requires  rn  inversions  of  I  —  p2 D/m,  and 
this  approach  can  be  interpreted  as  integration  of  the  DE  by  an  implicit 
scheme18  with  the  ’’time  step”  St  =  p2/m.  The  second  method  involves 
only  one  inversion  of  the  matrix  whose  condition  number  is  cm,  where 
c  =  cond(l  —  p2D/m). 
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Numerically,  each  iteration  of  the  first  inversion  process  is  approximately 
equivalent  to  a  one-step  integration  with  the  explicit  scheme  (7)  as  the 
corresponding  matrices  l-j-p2D/n  and  I— p2D/m  differ  only  by  the  numerical 
factor  before  the  diffusion  operator.  The  second  method  of  computing  the 
action  of  Bm  is  more  expensive  than  the  first  one,  because  the  total  number 
of  iterations  grows  exponentially  with  m,  unless  an  efficient  preconditioner 
is  available. 

On  the  other  hand,  the  possibility  to  directly  compute  the  action  of 
B-1  =  (I  -  D/m)”*  is  advantageous  in  solving  the  MS  3dVar  problem  (2), 
as  it  requires  only  one  MS  iterative  cycle  to  invert  (I  —  D/m)m  -f  H  1  H.  In 
contrast,  the  DS  solution  (4)  and  the  Bm-preconditioned  MS  solution  (5) 
involve  the  product  of  two  cycles:  each  iteration  of  the  respective  DS/MS 
system  solvers  contains  an  MS  iterative  cycle  required  for  computing  the 
action  of  B  (or  B)  on  a  vector. 

Spectral  properties  of  the  low-order  spline  models  differ  considerably 
from  the  Gaussian  one:  their  spectra  exhibit  more  gentle  slopes  and  weaker 
damping  of  the  short  (near-grid)  scales  (Fig.  la)  and  the  correlation  func¬ 
tions  decay  faster  than  the  Gaussian  at  small  distances  (Fig.  lb).  The  differ¬ 
ence  may  affect  the  forecast  skill  of  the  assimilation  system  and  not  worth 
the  computational  gain  when  applied  to  real  data.  This  and  other  related 
issues  have  been  examined  by  means  of  numerical  experimentation. 

3.  Numerical  experiments  setup 

Experiments  were  performed  with  the  Relocatable  Navy  Coastal  Ocean 
Model  system  (RNCOM)  consisting  of  two  primary  components:  The 
NCOM  provides  forecasts  of  the  ocean  state,  and  the  Navy  Coupled  Ocean 
Data  Assimilation  (NCODA)  uses  a  3dVar  algorithm  to  assimilate  obser¬ 
vations  into  the  model  forecast  state19. 


3.1.  Numerical  model  and  observations 

NCOM  has  a  free-surface  and  is  based  on  the  primitive  equations  under  the 
hydrostatic,  Boussinesq,  and  incompressible  approximations.  The  Mellor 
Yamada  Level  2/2.5  turbulence  models  are  used  to  parameterize  vertical 
mixing.  Most  terms  are  treated  explicitly  in  time,  except  for  the  propagation 
of  surface  waves  and  vertical  diffusion,  which  are  treated  implicitly.  For  the 
present  study  the  model  was  configured  on  two  grids  with  homogeneous 
grid  spacing  of  3  and  10  km  in  the  horizontal.  In  the  vertical  there  were 
respectively  46  and  50  layers  having  grid  steps  varying  between  1  m  and 
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Fig.  2.  The  model  domain  and  sea  surface  temperature  increments  at  10  m  on  Spetem- 
ber  2,  2007  for  the  C d  (a)  and  covariance  models  (b)  (see  Section  3.2).  Contour  interval 
is  0.1°. 


400  m.  The  number  of  grid  points  representing  a  3d  scalar  field  was  M  = 
10, 766, 576  and  M  =  862, 992  respectively.  All  the  runs  were  conducted  on 
the  Dell  R610  server  equipped  with  16  Xeon  5500  processors  running  at 
2.8GHz. 

Assimilation  experiments  were  performed  in  the  Okinawa  Trough  region 
(Fig.  2)  in  the  time  period  from  September  1  to  October  31  2007.  The  region 
and  time  period  were  selected  to  include  extensive  Navy  observations  from 
an  air-deployed  bathythermograph  survey,  a  shipboard  hydrographic  sur¬ 
vey,  and  eight  gliders.  Observations  from  this  Navy  exercise  are  an  addition 
to  the  standard  operational  data  stream  used  by  NCODA,  which  consists 
of  sea  surface  temperature  (SST)  and  sea  surface  height  anomalies  obtained 
from  satellites,  and  temperature/salinity  profiles  acquired  by  buoys,  floats, 
CTDs,  and  XBTs.  The  total  number  of  observations  processed  during  the 
2- month  assimilation  period  was  1,050,429,  or  approximately  N=  17,507 
points  per  24-hour  assimilation  cycle. 


3.2.  Assimilation  system 

NCODA  uses  a  DS  3dVar  data  assimilation  scheme  with  the  analysis  equa¬ 
tion  (4).  The  vector  of  analysis  variables  x  contains  temperature,  salinity, 
geopotential  (dynamic  height)  and  velocity  fields,  but  in  contrast  to  the 
DE  approach,  the  BEC  operator  is  defined  by  the  explicit  specification  of 
its  matrix  elements  via  correlations  in  3d  using  the  correlation  function 
shown  by  the  dashed  line  in  Fig.  lb.  In  the  following,  we  will  denote  this 
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BEC  model  by  C<*,  the  Gaussian  model  and  its  mth-order  spline  approx¬ 
imations  will  be  labeled  by  Coo  and  Cm,  and  the  asterisk  will  denote  MS 
implementation  (3)  of  the  spline  model  (C^). 

In  the  assimilation  experiments,  the  Cd  model  was  replaced  by  the  tested 
BEC  models.  Several  assimilation  runs  with  Cd  and  other  BEC  models 
were  also  performed  over  the  2-month  assimilation  period  for  comparison 
purposes.  In  these  runs,  the  horizontal  decorrelation  scale  was  set  to  45  km, 
while  vertical  scale  varied  in  z  in  proportion  with  the  vertical  model  grid 
step. 

Since  the  major  goal  of  the  present  study  is  to  compare  computational 
efficiencies  of  the  BEC  models  that  are  quite  different,  their  forecast  skill 
was  monitored  with  respect  to  the  operation  of  the  NCODA  system  with  the 
Cd  BEC  model  whose  forecast  skill  was  used  as  a  benchmark.  This  was  done 
to  ensure  that  the  computational  cost  of  the  analysis  was  not  reduced  at  the 
expense  of  reduction  in  assimilation  quality  q.  The  latter  was  estimated  as 
the  DS  distance  between  the  24-hour  model  temperature/salinity  forecast 
at  observation  points  7/,S/  and  the  observed  values  T0,S0: 

QT(t)  =  <(7>  -  T0)  V \T}  -  To))1'2,  (10) 

Qs(t)  =  ((Sj  -  S0)Tof(Sf  -  S0))1'2,  (11) 


where  <7t,s  are  the  observation  errors  and  the  angular  brackets  denote  aver¬ 
aging  over  the  observational  locations.  These  DS  distances  were  normalized 
to  measure  the  forecast  skill  s  of  the  tested  models  relative  to  the  skill  of 
the  benchmark  model  C<f. 


,.s  _  <]t(C)  +  qsjC) 
Qr(Cd)  +  Qs{Cd) 


(12) 


4.  Results 

4.1.  Comparison  of  the  forecast  skills 

As  it  has  been  noted  in  Section  2,  spline  models  are  characterized  by  broader 
spectra  and  provide  less  attenuation  at  high  spatial  frequencies  (Fig.  1) 
than  the  Gaussian  model.  This  property  causes  a  certain  difference  in  the 
analyses  increments  (Fig.  2),  which  may  result  in  substantial  decrease  of  the 
overall  forecast  skill.  The  forecast  skills  for  the  10  km  and  3  km  resolution 
configurations  are  shown  in  Figure  3.  It  is  seen  that  the  forecast  skill  of 
both  Coo  and  C2  BEC  models  does  not  depend  on  the  minor  changes  in  the 
shape  of  the  correlation  function:  The  2-month  mean  values  shown  in  Fig. 
3  do  not  differ  significantly  from  1.  This  result  indicates  that  the  analyses 
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Fig.  3.  Forecast  skill  of  the  C2  (black)  and  Coo  (gray)  BEC  models  implemented  on 
coarse  (10  km,  above)  and  fine  (3km,  below)  resolution  grids.  The  time-averaged  skills  s 
are  shown  in  the  low  right  corners. 


from  the  CQ 0  and  C2  models  are  at  least  not  degrading  the  benchmark 
forecast  generated  by  the  operational  NCODA  system,  while  both  models 
demonstrate  similar  forecast  skills. 

It  is  also  remarkable  that  the  forecast  skill  of  the  fine-resolution  models 
appeared  to  be  13-15%  below  the  skill  of  the  respective  coarse-resolution 
configurations.  To  some  extent  this  phenomenon  can  be  explained  by  the 
presence  of  small-scale  motions  in  the  3  km  configuration  that  are  barely 
constrained  by  the  available  observations:  On  average,  an  observation  sup¬ 
plies  information  for  610  grid  points  in  the  fine-resolution  case  against  110 
grid  points  per  observation  for  the  coarse-resolution  configuration. 


4.2.  Comparison  of  the  CPU  times 

The  dependence  of  CPU  time  was  explored  on  both  the  ratio  p  of  the  back¬ 
ground  decorrelation  scale  to  the  grid  step  and  on  the  degree  of  anisotropy 
of  the  correlations.  A  series  of  experiments  were  performed  with  differ¬ 
ent  strengths  of  the  anisotropy  and  different  values  of  p  for  the  selected 
date  September  2,  2007  (23,970  observation  points).  In  these  experiments, 
the  diffusion  tensor  was  specified  as  follows.  The  background  decorrelation 
scales  pi  at  every  location  were  defined  as  a  product  of  the  local  grid  steps 
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Fig.  4.  Composite  maps  of  four  columns  (filled  contours)  of  the  temperature  correlation 
operator  for  (a)  isotropic  (A  =  12)  and  (b)  anisotropic  (A  =  165)  cases  of  the  Coo  BEC 
model  at  20  m.  Thin  contours  show  the  SSH  field.  Contour  interval  is  10  cm. 


and  the  universal  scaling  parameter  p.  The  smaller  principal  axis  in  the 
horizontal  direction  (corresponding  to  A2)  was  set  to  be  orthogonal  to  the 
local  velocity  vector  v.  The  length  of  the  larger  axis  Ai  was  set  to  be  equal 
to  A2-max(l,  >J\v\/v),  where  v  is  a  prescribed  threshold  value  of  |v|.  A 
structure  like  this  simulates  enhanced  diffusive  transport  of  model  errors 
in  the  regions  of  strong  currents  on  the  background  of  isotropic  error  dif¬ 
fusion  (Fig.  4).  The  strength  of  anisotropy  was  controlled  by  changing  the 
value  of  v :  t/=10  m/s  corresponds  to  locally  isotropic  diffusion  (A=12,  Fig. 
4a),  v  =  0.2  m/s  imposed  moderately  anisotropic  covariances  (A  =  50)  in 
regions  of  strong  currents,  and  v  =  0.07  m/s  corresponds  to  the  strongly 
anisotropic  case  (A  =  165,  Fig.  4b). 

In  a  series  of  experiments,  NCODA  observations  on  Spetember  2,  2007 
were  analyzed  using  the  and  C  1,2,3  BEC  models  in  both  state-  and 
data-space  formulations  and  the  required  CPU  times  for  these  analyses 
were  compared.  Results  of  the  10  km  grid  size  experiments  are  assembled 
in  Table  1  where  larger  anisotropy  corresponds  to  the  larger  maximum 
eigenvalues  A  of  the  diffusion  operator  (column  2).  Respectively,  the  tested 
values  of  p  correspond  to  the  decorrelation  scales  of  30,  45  and  70  km. 

As  can  be  seen  from  Table  1,  the  numbers  indicate  improved  computa¬ 
tional  efficiency  of  the  low-order  (m  <  3)  MS  implementation  of  the  spline 
model.  For  higher-order  models,  the  MS  solution  appears  to  be  less  efficient 
due  to  exponential  growth  of  the  condition  number  of  the  system  matrix. 
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Table  1.  CPU  times  Too  for  the  Gaussian  covariance  model  and  the  relative 
CPU  times  of  the  mth-order  spline  models  implemented  in  DS  (rm)  and  in 
MS  (r,^).  The  accuracy  of  system  solutions  (defined  as  the  ratio  between  the 
norms  of  the  residual  and  the  rhs  vectors)  is  e  =  10'“6.  The  fastest  cases  for 
a  given  m  are  boldfaced. 


p 

A 

Too,  min 

rm  /  Too 

Tm/Too 

Coo 

m=  1 

m=  2 

m=3 

m=l 

m— 2 

m= 3 

12 

3.47 

0.83 

2.00 

2.66 

0.06 

0.31 

1.36 

3.0 

50 

14.9 

0.38 

0.85 

1.10 

0.O2 

0.21 

1.82 

165 

54.0 

0.16 

0.32 

0.40 

0.01 

0.15 

1.78 

12 

8.72 

0.54 

1.25 

1.54 

0.04 

0.23 

1.58 

4.5 

50 

36.4 

0.19 

0.41 

0.74 

0.01 

0.15 

1.62 

165 

156 

0.09 

0.21 

0.29 

0.00 

0.09 

1.83 

12 

24.0 

0.24 

0.62 

1.06 

0.02 

0.16  1 

1.42 

7.0 

50 

129 

0.08 

0.19 

0.27 

0.00 

0.07 

0.84 

165 

418 

0.04 

0.11 

0.16 

0.00 

0.06 

6.68 

Additional  experiments  were  made  on  a  longer  time  scale,  with  the 
NCODA  system  using  the  generic  correlation  Cd  model  as  well  as  the  tested 
models  with  24-hour  analysis  cycle  (A  =  12,  p  =  4.5).  These  experiments 
have  shown  that  the  CJ  model  is  3  times  faster  than  C2  and  3.5  times  faster 
than  Co©  for  the  10  km  configuration.  Similarly,  for  the  3  km  configuration, 
the  C 2  model  was  3.3/4. 2  times  faster.  CPU  times  of  the  C £  and  Cd  models 
are  compared  in  Fig.  5.  On  average,  the  C £  model  requires  30-50%  more 
CPU  time  than  the  generic  Cd  model.  However,  when  the  number  of  ob¬ 
servations  exceeds  1.5-1.71 04 ,  the  C £  model  appears  to  be  more  efficient. 
Similar  computations  for  3  km  resolution  show  that  this  critical  number  of 
observations  increases  only  slightly  to  1.8-2*  104  despite  a  12-fold  increase 
in  the  dimension  of  the  model  space. 


5.  Conclusions 

The  forecast  skill  and  computational  efficiency  of  the  Gaussian  and  spline 
covariance  models  were  examined  in  the  framework  of  3dVar  assimilation 
of  real  data  into  an  operational  ocean  model.  It  is  shown  that  the  MS  for¬ 
mulation  of  the  second-order  spline  model  has  similar  24-hour  forecast  skill 
and  3-5  times  better  computational  efficiency  than  the  DS  implementation 
of  the  Gaussian  and  spline  models.  At  m  <  3,  the  computational  efficiency 
of  the  C m  solutions  is  based  on  the  low-cost  computation  of  the  action  of 
the  inverse  BEC  operator  B  1  =  (I  —  p2D/m)m  which  contributes  to  the 
system  matrix  of  the  normal  equation  (2).  On  the  contrary,  multiplication 
by  I  —  p2 D/m  (or  by  I  *f  p2D/n)  has  to  be  performed  many  times  in  the  DS 
formulation  to  iteratively  model  the  action  of  B,  which  in  turn  is  immersed 
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Fig.  5.  CPU  time  ratios  between  the  Cj  and  C\  models  against  the  number  of  obser¬ 
vations  for  the  10  km  configuration.  Points  above  the  horizontal  line  correspond  to  the 
cases  of  the  faster  model. 


into  the  iterative  loop  required  to  find  the  solution  (4)  to  the  normal  system 
in  the  DS  formulation. 

It  is  also  shown,  that  the  difference  in  the  BEC  models  has  a  negligible 
impact  on  the  forecast  skill  of  the  3dVar  assimilation  system.  Comparison 
with  the  benchmark  NCODA  3dVar  algorithm  has  shown  that  the  forecast 
skill  remains  virtually  the  same  (Fig.  3),  whereas  the  model  appears  to 
be  more  efficient  computationally  than  the  operational  BEC  model  when 
the  number  of  observations  exceeds  15- 12- 103.  Numerical  experiments  have 
also  shown  that  spline  models  become  especially  advantageous  when  the 
background  decorrelation  scale  is  well  resolved  by  the  model  grid  (p  >  3) 
and  the  diffusion  operator  is  strongly  anisotropic/inhomogeneous  (Table  1). 

The  results  of  this  work  suggest  that  studying  the  applicability  of  the 
anisotropic  higher-order  spline  BEC  models  to  3dVar  assimilation  is  worth 
consideration  for  at  least  three  reasons:  1)  they  are  computationally  efficient 
in  processing  large  number  of  observations,  2)  they  are  flexible  enough  to 
accommodate  covariance  information  from  the  structure  of  the  background 
flow,  and  3)  they  can  be  easily  extended  to  include  model-generated  covari¬ 
ances  extracted  from  model  statistics. 
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